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Abstract 

The Four Fermi model with discrete chiral symmetry is studied in three dimensions at non- 
zero chemical potential and temperature using the Hybrid Monte Carlo algorithm. The number 
of fermion flavors is chosen large {Nf = 12) to compare with analytic results. A first order chiral 
symmetry restoring transition is found at zero temperature with a critical chemical potential ^c in 
good agreement with the large Nf calculations. The critical index v of the correlation length is 
measured in good agreement with analytic calculations. The two dimensional phase diagram (chem- 
ical potential vs. temperature) is mapped out quantitatively. Finite size effects on relatively small 
lattices and non-zero fermion mass effects are seen to smooth out the chiral transition dramatically. 
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1. Introduction 

The behavior of symmetries at finite temperatures and densities is one of the most outstanding and 
relevant problems in many current areas of particle physics; e.g. cosmology, rclativistic heavy-ion collisions, 
and the quark-gluon plasma [1] . In recent years we have witnessed revived interest in the chiral symmetry 
restoration transition in QCD. The problem of symmetry breaking and its restoration is intrinsically non- 
perturbative. Therefore, the number of available techniques is limited and most of our knowledge about the 
phenomenon comes from lattice simulations. Due to the enormous complexity of QCD, studies have so far 
been done on lattices of modest size and have been unable to yield quantitative claims as far as the order 
of the transition is concerned. This is unfortunate since several studies suggest that the high temperature 
phase of QCD has a number of interesting features [2]. In addition, only very slow progress has been made 
in lattice simulations at finite chemical potential, which is the regime of direct relevance to real physics [3] . 

In this paper we will approach the general problem of chiral symmetry restoration at finite temperature 
and density in a three-dimensional toy model in order to understand what ingredients might play a decisive 
role in more complex systems like gauge theories. We have simplified our model as much as possible in order 
to produce the highest quality data and learn what range of parameters we need for studies of more realistic 
cases. We have thus chosen to study the Gross- Neveu model [4], in which the chiral symmetry is discrete. 
The Lagrangian is 



Nf 



^0)^^0)_ J^(^0)^0))2 



(1.1) 



Here '0 is a four component spinor, Nf is the number of flavors of elementary fermion, and the discrete 
Z2 chiral symmetry is "0 '"*' TsV'i '~* "V'Ts- Our choice of model also reflects our interest in its behavior at 
zero temperature and density [5,6]: in less than four dimensions the model has a non-trivial renormalization 
group fixed point, also characterized by the spontaneous breaking of the chiral symmetry at a critical coupling 
g^. It is thus a toy model for the study of non-trivial strongly coupled theories. Although the theory is 
non-renormalizable in a standard perturbation expansion in the coupling g^, its 1/Nf expansion about the 
fixed point gj is renormalizable. We have argued elsewhere [6] that this is closely related to the fact that 
the theory's critical indices satisfy hyperscaling. Physically this means that the theory has a divergent 
correlation length at the fixed point and this length sets the scale for the theory's low-energy phenomena. 
To 0{1/Nf) the theory's critical indices are [6] 

''='^ + T^' ^ = ^ + Tr^'^ /^'"^l; 7 = l + T7^; ^ = l-:71^- (1-2) 

in the standard notation of classical statistical mechanics [7] . 
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There are several motivations for studying such a simple model. In four dimensions the four-fermi model 
is believed to be an effective theory of quarks and gluons at intermediate energies. The degrees of freedom 
are light mesons and quarks. As far as finite temperature and density is concerned, the low temperature 
regime will be dominated by the lightest particles and, if the restoration temperature is of the order of 
lOOMeV, then the contribution of the heavier particles like p mesons will be exponentially suppressed. In 
that sense, the universal properties of chiral symmetry restoration in QCD could be well described by an 
effective theory like the Nambu ~ Jona-Lasinio model. In three dimensions, however, the four-fermi model is 
actually renormalizable, and precise analytic predictions are available from an expansion in l/Nf. Although 
we consider the simplest such model in this paper, extensions to more realistic models where the chiral 
symmetry is continuous are straightforward [4]. One might think of a yet more drastic simplification, and 
consider the model in two dimensions. However, in this case there are conceptual problems; eg. in the Z2 
case the symmetry restoration is now dominated by the materialization of kink - anti-kink states, which 
are composite states of the fundamental fermion fields, which arc not probed in the l/Nf expansion. In 
two dimensions the extension to continuous symmetries is also plagued by the non-existence of massless 
Goldstone bosons [8]. 

The major technical barrier to progress in simulating non- vanishing baryon number densities in QCD is 
the absence of a probabilistic interpretation of the path integral measure due to the action becoming complex 
once the chemical potential fi ^ 0. In the model considered here, the action remains real even after the 
introduction of /i, which means we can study the physics of the high-density regime using standard Monte 
Carlo techniques. In addition, it has fewer degrees of freedom than QCD, or even the Nambu - Jona-Lasinio 
model, and hence can be studied with greater precision on much bigger lattices than are presently used for 
QCD thermodynamics. In general, simulations focus on the evaluation of the order parameter < ipip > which 
requires the inversion of the Dirac operator. Most numerical algorithms which simulate the effect of virtual 
quark anti-quark pairs in the vacuum also require the inversion of this operator; this is the most computer- 
intensive step in such simulations. In gauge theories, this operator is singular in the chiral limit (ie. the 
matrix to be inverted becomes ill-conditioned), and the simulations have to be done using a finite bare mass. 
Information about the chiral limit is then obtained by extrapolating the finite mass data to ?Ti = 0. This last 
step is severely constrained by the lattice volume; m can not be taken arbitrarily small, since otherwise the 
Compton wavelength of the Goldstone pion associated with the symmetry breaking would exceed the lattice 
size producing severe finite volume effects [9]. Consequently, lattice QCD data, always taken at finite mass, 
show only a crossover rather than a real phase transition. In our model, we are dealing with a Yukawa-like 
coupling (see below), so that the fermion matrix to be inverted has non-zero diagonal elements even when 
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the bare mass is set to zero. The order parameter (the inverse Dirac operator) is thus not singular and the 
simulations can be done in the chiral limit directly. In this way, we can in principle explore the systematics 
of finite mass effects, finite volume effects, and the validity of various extrapolation procedures. This might 
teach us how to use the QCD data better. 

The mean field theory description of the transition, to be presented in Sec. 3, predicts a first order 
transition for T — and a continuous transition for T > 0. Both analytic and numerical work in this model 
is aided by the introduction of an auxiliary scalar field a, so eqn. (1.1) becomes 

Nf 



C^J2 fV'^^'-'^V^^^' + crVS^^'V^^'^1 + T-{a\ (1.3) 

and the Lagrangian becomes quadratic in the fermion field. The ground state expectation value of a serves 
as a convenient order parameter for the theory's critical point. Using the standard lattice regularization 
scheme to be discussed in Sec. 2, the bulk critical coupling is found to be 1/g^ — 0.975 - 1.000 for Nf — 12 
when the chemical potential is set to zero as in eqn. (1.3) [5]. Simulation of 20'^ lattices at various couplings 
close to the critical point (within its scaling region as determined in ref. [5]) showed that S =< a > jumps 
discontinuously to zero as /i is increased and that the induced fermion density also jumps dis continuously 
through the transition. In addition, the dependence of ^c on the coupling allows us to determine the 
correlation length exponent v for Nf = 12: 

1^ = 1.05(10), (1.4) 

in good agreement with the large Nf prediction (1.2). Furthermore the magnitude of ^c itself is found to 
be in good agreement with the mean field result /ic = So , where Eg is the value of the vacuum expectation 
of the scalar field at zero temperature and chemical potential. This result indicates that materialization of 
the fermion itself drives the symmetry restoration transition: this is not the case in the two-dimensional 
Gross- Neveu model, where kink - anti-kink states materialize at the transition [10,11]. The fact that ^c 
agrees with the mean field result is indirect evidence that such exotic fermion states in which the energy per 
constituent is smaller than the energy of a single fermion state do not occur in the three dimensional model 
[12]. 

We also simulated the model at both non-zero /i and non-zero temperature T. In ref. [5] we confirmed 
that the /^ = 0, T 7^ symmetry restoring transition occurs at Tc/Tio ~ 0.72 in good agreement with mean 
field predictions, and that this transition is second order. Here we map out the phase diagram in the (/i, T) 
plane. 

The simplicity of this model and the efficiency of the simulation algorithm allowed us to address two 
technical issues of interest to lattice gauge theorists studying QCD in extreme environments. First, four- 
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dimensional QCD simulations are presently restricted to relatively small lattices. We show here that if the 
lattice is taken relatively small (12"^ as opposed to 20^) the non-zero /i transition is smeared out and all 
evidence for a discontinuous transition is lost. Furthermore, as discussed above, QCD cannot be simulated 
directly in the chiral limit. We show here that, even on relatively large lattices and at couplings chosen in 
the scaling region, even very small bare fcrmion masses m smooth the transition dramatically. We learn that 
it would be very difficult indeed to find evidence for a discontinuous transition by approaching the chiral 
limit TO = of the theory via to. 7^ simulations. 

This paper is organized into several sections. In Sec. 2 we briefly discuss the lattice formulation of the 
model and its simulation algorithm. More detail on these issues has already been provided in ref [5]. In Sec. 
3 we present the mean field analysis of the theory - this overlaps to some extent with the results of ref. [12]. 
In Sec. 4 the simulation study of the zcro-tcmpcrature model at non-zero chemical potential is presented. In 
Sec. 5 the model is sudied at both /i 7^ and T 7^ 0. Finally, in Sec. 6 we consider finite volume effects by 
simulating on a 12'^ lattice, and finite m effects on relatively large lattices. In both cases the discontinuous 
nature of the transition driven by non-zero chemical potential is lost. Sec. 7 summarizes our work. 

2. Lattice Formulation of the Gross-Neveu Model 

The Gross-Neveu model in its bosonised form (1.3) may be formulated on a space-time lattice using the 
following action: 



^//2 / ^ \ N- 

i—l y x.y X <x,x> j X 



(2.1) 



where XiiXi ^re Grassmann- valued staggered fermion fields defined on the lattice sites, the auxiliary scalar 
field a is defined on the dual lattice sites, and the symbol <x^x> denotes the set of 8 dual lattice sites x 
surrounding the direct lattice site x. The lattice spacing a has been set to one for convenience. The fermion 
kinetic operator M. is given by: 



J^x.,y = 2 






-^^ T]y{x){8y^xJ^0-iy,x~v\, (2.2) 



2 

1^=1,2 



where ri^{x) are the Kawamoto-Smit phases (— 1)^"^ \-x,,-i_ ^j-^g influence of the chemical potential fi is 
manifested through the timelike links, following [13]: only fermion loops which wrap around the timelike 
direction are affected by its inclusion. 

The Gross-Neveu model in two dimensions was first formulated using auxiliary fields on the dual sites 
in reference [14]. We can motivate this particular scheme by considering a unitary transformation to fields 
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u and d [15]: 

^ (2.3) 

Here Y denotes a site on a lattice of twice the spacing of the original, and ^4 is a lattice vector with entries 
either or 1, which ranges over the corners of the elementary cube associated with y, so that each site on 
the original lattice corresponds to a unique choice of A and Y . The 2x2 matrices V ^ and B^ are defined 

by 

(2.4) 



where the Th are the Pauli matrices. Now, if we write 

and interpret g as a 4-spinor with two flavors counted by the latin index a, then the fcrmion kinetic term of 
the action (2.1) may be recast in Fourier space as follows: 

——3^ ^ -<^ gi(fc)(7^(g)l2)gi(fc)sin2fc^ + gj(fc)(75 «)T*)qi(fc)(l-cos2fci.; 



i i/=l,2 



--< gi(fc)(7o ® \2)qi{k) [isin2/co cosh/i + (1 + cos2fco) sinh/i] (2-6) 

+9i(fc)(75 ® ^3 )9i(fc) [^(1 ~ cos2fco) cosh^ + sin2fco sinh/i] 



where 

{i.)c.p ^(^'' ) ; (70).^ ^(^' ) ; h,U = f ,1 ~'^0 , (2.7) 

the second set of (2 x 2) matrices in the direct product act on the flavor indices, and the momentum integral 
extends over the range fc^ G (— 7r/2,7r/2]. At non-zero temperature the lattice has finite extent in the 
temporal direction, and J dko is replaced by a sum over Nr/2 modes, where Nr is the number of lattice 
spacings in the time direction, and antiperiodic boundary conditions are imposed on the fcrmion fields. In 
the classical continuum limit lattice spacing a ^ 0, the fiavor non-diagonal terms vanish as 0(a), and we 
recover the standard Euclidian form qj{<^ + lJ-ja)Qj, where the flavor index j now runs from 1 to Nf. 

Similarly, it is straightforward to show that the interaction terms can be rewritten (with obvious nota- 
tion): 

" " ~ q^{Y){U®l2)q^{Y) + 0{a), (2.8) 



^»*-e(E'^(^;^)) 



where the 0{a) terms contain non-covariant and flavor non-singlet terms. If we used a formulation in which 
the a fields lived on the direct lattice sites, then such non-covariant terms would contribute at 0(a") [14]. 
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Thus we see that the lattice action (2.1) reproduces the bosonised Gross-Neveu model at non-zero 
density, at least in the classical continuum limit. Most importantly, (2.1) has a discrete global invariance 
under 

X^ix) ^ (-l)"«+"^+"^X^(^); X^(^) ^ -(-l)""+"^+"^X^(a:); a(x) ^ -a(x), (2.9a) 

ie. 

q^{Y) ^ i-r^ ® l)q,{Yy, g,(y)^-g,(r)(75(»l); a(£)^-a(x). (2.96) 

It is this symmetry, corresponding to the continuum form ip h^ 75'0, ip '~^ ^"075 7 which is spontaneously 
broken at strong coupling, signalled by the appearance of a non-vanishing condensate < XX > or equivalently 
< q{l4 <E) 12)? >• We shall see in the next section that to leading order in l/Nf the lattice formulation (2.1) 
gives predictions in agreement with those of the continuum. 

The action (2.1) was numerically simulated using the hybrid Monte Carlo algorithm [16], in which the 
Grassmann fields are replaced by real bosonic pseudofermion fields (?!)(x) governed by the action 

^ = E E 2<^^(-)(^'^):,V.(y) + i3E-'(^)' (2-10) 

x.y i,j=l ^ X 



where 



Mxyij ^ MxySij +5xy5ij- ^ (7{x). (2-11) 



<a:,a:> 



Note M is strictly real. Integration over <}> yields the functional measure A/det(M*M) = detM if the 
determinant of M is positive semi-definite. This condition is fulfilled if Nf/2 is an even number, even for 
/i 7^ 0. The problem of complex determinants associated with simulating gauge theories at finite density do 
not appear. Further details of our simulation procedure are given in [5]. 

As well as measuring the expectation value of the scalar field < cr > in the simulation, which for our 
purposes is the order parameter of the transition, we also monitored the chiral condensate < ipip >, the 
energy density < e >, and the fermion number density < n >, which are defined by 

- < ■0V' > == Y^'^^F = 77 < trM-i >, 

1 ainZ 1 „ 1 V- „ 1 „ 1 

< e >= -tf^T:^ = T7tr<9o7oS'F = TTTT < > e^M^^ - - e^^M^^ - >, .„-,„^ 

V, 9/3 V 2V ^-^ x.x+O a;,a;-0 (.2-12J 

X 

VsP da V 2V ^-^ a:,a:+0 x,x-l3 

X 

Here Vg is the spatial volume, /3 the inverse temperature, and V = Vsf3 the overall volume of spacetime. The 
final expression in each case is the quantity measured in the simulation, using a noisy estimator to calculate 
the matrix inverses. 
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3. Mean Field Analysis 

In this section we calculate the order parameter E =< cr > as a function of coupling 5, chemical potential 
/i, and temperature T (= 1//3), to leading order in 1/Nf. Since the limit Nf -^ 00 suppresses fluctuations 
around the saddle point solution, this is equivalent to a mean field treatment. We will work both in the 
continuum and using the specific lattice regularisation (2.1). If diagrams of 0{l/Nf) and beyond are ignored, 
the only contribution to S comes from a simple fermion loop tadpole, and we determine S self-consistently 
using the gap equation: 

j: = -g^ <{^^b>^^tTSF{fi,T,^), (3.1) 

where in the Euclidean formulation the fermi propagator Sp is given by 

Sp^{k;n,T,Y.)=ijo{ko-iiJ.)+ ^ ifc^7^ + S. (3.2) 

For non-zero temperatures the alUowed values of fco are quantised as 

ko ^ (2n - l)7rT, neZ, (3.3) 

ie. with antiperiodic boundary conditions in the finite temporal direction. Collecting together equations 
(3.1-3) we arrive at 

l = 4rV/^ '—2 • (3-4) 

5' „f"oo^ (2^)2((2n-l)7rr-ZA*)'+p2 + S2 

The manipulations from here are standard [10,12]. First one resums over n using the Poisson formula: 

1 V^ /" d'^P /"°° g27rim0 



rn— — 00 



Here E = ^jp- + Y? . If /i < S, then the poles in the integrand lie on opposite sides of the integration contour 
for all values of p, and the integral over </) is easily performed to yield 
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1 ("^ dE 

"2 



1- ^ 



(3.6) 



If, however, /i > S we must take care, since for certain values of p, both poles in the integrand of (3.5) will 
lie to the same side of the contour. We find in this case 



1 f°° dE 



9^ Jf, T^ 



1 
1 



gPiE-i,) _^_ I 



f^dE 1 rdE 1 ^g^^ 



/s TT e'3(^--^) + l 7s 71" e/^^^+^' + l' 
We now eliminate g in favour of So, the order parameter at zero temperature and chemical potential, using 
equation (3.6) at T = /i = 0. This gives an implicit equation for S in terms of a physical scale Eq; with no 
reference to any UV cutoff. For /z < S we obtain 

So - S = T (ln(l + e-'^f^-^)) + ln(l + e~^'^^+^'^)^ . (3.8a) 
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while for /i > E: 

So - A* = r (ln(l + e-'3(A'-s)) + 1^(1 + g-/3(p+s))^ ^ (3 g^^ 

In fact, these two equations are identical solutions for S(/i, T), and also demonstrate that curves of S(/i) 
at constant T are symmetric under reflection in the line S = /i. This result is peculiar to three spacetime 
dimensions. Equation (3.8a) was first derived in [12]. 

Equation (3.8) gives a complete solution for Y,{ji,T) in terms of Sq. Since E ^ 0+ smoothly, the 
symmetry-restoring transition is second order throughout the (/i, T) plane, except for one isolated point, as 
we shall see. To obtain the equation for the critical line in this plane we set S = in (3.8) to get the curve: 

l-^ = 2^Hl + e-^n. (3.9) 

At zero chemical potential, therefore, we predict a chiral symmetry-restoring transition at a critical temper- 
ature 

T. = ^ ^ 0.72E„. (3.10) 

The gap equation in the broken phase in this limit is the /i = limit of (3.8): 

Eo-E = 2Tln(l + e"''^). (3.11) 

At zero temperature we find E = Eg independent of /i up to a critical value 

^^c = Eo, (3.12) 

at which point there is a discontinuous drop to zero, ie. at this isolated point the transition is first order. 
For small excursions into the (/i, T) plane we find from (3.8) 

|^|T^o=nm-^.-e^(--), (3.13) 

d^i T-,0 smh/3E 

ie. the slope of the surface E(/x,T) diverges in an essentially singular way as ^ -^ ^c, T -^ 0. So, the mean 
field analysis predicts a first order transition for T = 0, which becomes second order as soon as T > 0. 

Using a similar route we can also calculate the fermion number density < n > in the broken phase 
starting from (2.12). For /i < E we find 

ET, (l + e-^(^"'')) 2T2^, ,.,e-^'^^sinh/3fcAi . ^ 

TT (l-he-'5(^+^)) TT ^^ ' k^ ^ ' 

k—1 

whereas for /i > E: 

^2-S2 ET, (l + e-^(^-^)) 2T2^, , Jl-e-^'^^cosh/3/cE) 



< n >= 



'° :-«.. -^E(-')' ^"' .."'"' ■ ("5) 



2tT TT (1+6-/5(^+5^)) TT ^^^ ^ fc2 

^ ^ k=l 
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In the symmetric phase we recover the usual expression for a two-dimensional relativistic free fermi gas: 

In the limit T ^ we see that ferniion density is strongly suppressed for /i < S, then jumps discontinuously 
and continues to rise quadratically with /x as soon as ^ exceeds E. As required, < n > vanishes for all T 
when 11 — 0. 

We have also studied the gap equation using an explicit UV regularisation defined by the lattice action 
(2.1). This is useful for comparison with the numerical results, particularly so that we can determine whether 
any differences with the continuum predictions arise from genuine l/Nf corrections (ie. departure from mean 
field behaviour) , or simply from the fact that on a finite lattice it is impossible to attain the thermodynamic 
limit. This latter point arises because lattice simulations at non-zero temperature are generally accomplished 
using a system with Nr lattice points in the temporal direction, with A^,- ^ N , the spatial dimension. In 
the work presented here N^- ranges from 6 to 12, for iV = 36. Clearly the main effect is that the sum over 
Matsubara frequencies in the lattice gap equation is truncated at a rather small value of n. 

Using the free fermion action in the form (2.6) we evaluate equation (3.1) to yield the lattice gap equation 
on a system of infinite spatial extent: 

8 r/^ (Pk ^ 1 






cos 



a . /27r(2n- 1, , . , ^ 

2'^H^ — '^^^ ^ 



■^ 1 — cos I -^-^ — - ) cosh 2/i — 2 sin ( —^ — - ) sinh 2/i > + '}Zv=i 2 ^^^^ ^"^ + ^^ 

(3.17) 

where ^^ defines a sum running from — ^ + ito-^ — ^if Nr/2 is odd or — ^ + 1 to ^ if iV^/2 is even 

(note Nt must be even in order for staggered fermions to be defined). The integral over k can be done via 

Schwinger parameterisation to yield 

1 = ^/°°c^ae-"^i+-^)/o^(f)E-P [f cos (^^^) cosh2;." 

(3.18) 
which is now in a form suitable for numerical quadrature, /o is the modfied Bessel function. Unfortunately 
the integral over a in (3.18) only converges for 2/i < cosh~ (sec27r/A^T-): for values of /i greater than this, 
although (3.17) is convergent, it is no longer possible to cast it into convenient form. 

Figure 3.1 shows a comparison of T,{p) calculated using the continuum solution (3.8) and the lattice 
solution (3.18) evaluated at inverse coupling 1/g^ = 0.75 on a lattice with Nr — 6, 8, 10 and 12. In practice we 
solve (3.18) for 1/ g^ as a function of S and invert using interpolation. The continuum solution sets (3 — Nr 
and uses a value for Sg given by the lattice gap equation at zero temperature and chemical potential: 

9 Jo 
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We see that the agreement between the two is fair, with the lattice rcsuhs always lying below the 
continuum ones. As discussed above, we ascribe this difference to the finite number of thermal modes 
available on the lattice. 

4. Non-zero Chemical Potential Simulation Results 

We first studied the theory at non-zero chemical potential on a symmetric lattice. In ref. [5] we 
obtained accurate results at vanishing chemical potential for lattices ranging in size from 8'^ through 20"^. 
The vacuum expectation value of a, and its susceptibility were measured over a range of coupling extending 
from 1/g^ ^ 0.5 to 1.2. Good agreement with large Nf scaling laws were found for 1/g^ ranging from 0.70 to 
1.1 and the critical point in the infinite volume limit was estimated to be l/g^ — 0.975—1.000. So, the scaling 
window of the lattice formulation lying in the chirally asymmetric phase was seen to be 0.70 < 1/g^ < 1.00. 
These results led us to simulate the model with /i 7^ on a 20^ lattice at couplings 1/g^ = 0.70,0.75 and 
0.80. Short exploratory runs indicated that larger lattices would be needed to push l/g"^ closer to the critical 
point. In Table 1 we show the simulation results for various /i at l/.g^ = 0.70, measurements of S =< a >, 
the energy density < e >, and the induced ground state fermion number density < n > are recorded. We 
also show the number of trajectories of the hybrid Monte Carlo algorithm at each point. The huge number 
of trajectories relative to state-of-the-art lattice QCD simulations with dynamical fcrmions was possible 
because of the three dimensional character of the model and its relatively simple form, a random scalar field 
coupled to a fermionic scalar density. This last feature led to a conjugate gradient routine which converged 
with an order of magnitude fewer sweeps than typically needed in lattice QCD simulations. The results 
recorded in Table 1 arc plotted in Figure 2 where we see a jump discontinuity in S vs. /i of 0.275 as /i varies 
from 0.39375 to 0.3941. Note from the table that we were particularly careful to accumulate good statistics 
near the transition. In Figure 3 we show the induced fermion number < n > plotted against /i which shows 
a clear jump over the same coupling range. We believe that the nonvanishing values of < n > recorded 
for fi < 0.3941 are finite size (temperature) effects. Finite size effects will be discussed further in Section 
6 below. The error bars in the table and plotted in the figures come from standard binning procedures. It 
was possible to "confirm" these error estimates in many cases by running the algorithm for an extra several 
thousand trajectories and reproducing average values and variances. 

As shown in Table 2 the simulation was repeated at 1/g^ = 0.75, slightly closer to the bulk critical 
point. The critical chemical potential shifted to iic ~ 0.32 — 0.3225 and discontinuities were again observed 
in S vs. /i (Figure 4) and < n > vs. /i (Figure 5). The discontinuities were slightly more difficult to measure 
because of the proximity to the critical point which reduces the size of physical observables when measured 
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in lattice units. 

Runs were also completed at 1/g^ ~ 0.80 and a value of ^c. ~ 0.250(5) was measured. Less extensive 
data was taken here as reflected in the larger error bar. Wc identified [i^, by plotting time histories of S (S 
vs. computer time) and noting that for ^ < jic maintained a non-zero value while for fi > iic evolved to zero 
and fluctuated around it. These results were less quantitative than the l/g^ = 0.70 and 0.75 simulations 
because critical slowing down was affecting the efficiency of the algorithm. In fact, we abandoned an attempt 
to study l/g^ — 0.85 because of severe critical slowing down and the possibility of large finite size effects 
invalidating an estimate of fic on a 20'^ lattice. 

It is interesting to analyze the measurements of ^c at l/g"^ = 0.70, 0.75 and 0.80 for two additional 
purposes. The first is to estimate fic in physical rather than lattice units. The natural way to do this is to 
record the ratio of fic to Sq measured at the same coupling as the finite-m transition. At large Nf, Eq is 
essentially the dynamical fcrmion mass and Mean Field theory predicts that /ic/^o = 1.0. From rcf. [5] we 
have the values of Sq on a 20^ lattice (^ ^ 0) at l/g^ = 0.70(i;o = 0.432), l/g'^ ^ 0.75(1)0 = 0.346) and 
1/g^ ~ 0.80(So — 0.262). The ratios /ic/So at each coupling are plotted in Figure 6. The error bars come 
almost exclusively from the ^c measurements - the 20"^ /i = measurements of Sq recorded in rcf. [5] were 
very accurate indeed. The l/g^ = 0.80 result for ^c is particularly uncertain. Nonetheless, as the critical 
point is approached the curve strongly suggests that /ic/^o approaches unity in accord with Mean Field 
theory, although a departure downwards from this value as a result of l/Nf corrections cannot be ruled out. 

Our last use for this data is for a calculation of the critical index v, the critical index of the correlation 
length. Since the critical chemical potential fic is a dimensionful parameter coupled to a conserved current 
in the Lagrangian, it undergoes no rcnormalization due to l/Nf corrections, and should scale as a physical 
mass as the critical point is approached, i.e. vanishing in lattice units with the exponent i^ 

l^c = C{llgl-llg^y (4.1) 

From rcf. [5] I/5J = 0.975 — 1.00. In Figure 7 we show a plot of In /ic vs. ln(l/g^ — l/g1) with choice 
1/.9* — 0.975. A linear fit is good and it gives v ~ 1.00(5) to be compared with the first two terms of 
the large-TV/ expansion {v = I + ^/'iNfir'^ = 1.0225.., for Nf = 12). A similar plot for l/gl = 1.00 gives 
v — 1.10(5). So, in summary wc have 

J/ =1.05(10) (4.2) 

in excellent agreement with the \a,Tge-Nf analysis. This analysis and Eq. (4.2) could certainly be pursued 
more systematically with greater control, but consistency between the analytic and numerical approaches to 
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this problem is our only goal here. 

A systematic search for \/Nf effects in our simulation results, or in more accurate simulations on larger 
lattices that might be done in the future, would require more calculations than done here. For example, if 
we take our results for /ic/^o at face value, ie. an increasing function of l/.g^ as l/.g^ -^ l/5*i then we would 
predict the exponent v to be smaller than the exponent Pm governing the scaling of Sg, independently of 
any estimate of the exact value of l/gl- This is at variance with the 0{\/Nf) corrections of eq.(1.2) which 
predict 

i^M _ {i/gl-i/g^r _ ,^, 2 ^i„2.in^ .4 3^ 

So(5) Hlgl^ilg'Y- ^"* "^ ^ ^ 

Numerically 8/3NfTr'^ — 0.0225 for Nf == 12, and since 1/g^ — 1/g^ varies from 0.3 to 0.2 over the region of 
couplings explored here, the right hand side of eq.(4.3) varies by less than a percent. Clearly much greater 
precision is needed before we can interpret the trend of Figure 7 to l/Nf effects, which would require large 
corrections at Oil/N'j) and beyond. 

5. Non-Zero Chemical Potential and Temperature 

In ref.[5] we studied the four Fermi model at non-zero temperature by simulating it on asymmetric 
lattices Nr x iV^ with A^ 3> Nt and Nr ranging from 4 to 12. A second order chiral transition was discovered 
with a critical temperature measured in physical units of Tc/Sq ~ 0.70(5). The result was in good agreement 
with Mean Field theory which predicts Tc/T,o = (2 In 2)^^ = 0.721. This result and the zero temperature, 
non-zero chemical potential result of Section 4, were then extended to map out the /i — T phase diagram. 
We simulated 6 x 36^, 8 x 36^, 10 x 36^ and 12 x 36^ lattices and obtained mc on each of these lattices by 
measuring S vs. /i. Each lattice was simulated at a coupling 1/g^ = 0.75 because this coupling lay in the 
scaling window of the lattice Lagrangian and because of the success of Section 4. The curves of S vs. /i are 
shown in Figure 8. The statistics at each point are comparable with Table 1 and the error bars on each point 
in the figure are smaller than the symbols themselves. The precise values of fXc for each A^^- are recorded in 
Table 3 with error bars. As usual, vacuum tunnelling (flipflops between ±1] values) limited our precision 
on mc measurements. It is interesting to compare Fig. 8 with its Mean Field counterpart in Fig. 1. The 
qualitative features of both plots are identical, but the scale of S and fi values coming from the computer 
experiment are consistently below the Mean Field predictions. 

The results in Table 3 can be converted to physical points on a /ic/5^o vs. Tc/Sq phase diagram by 
recalling the value of So at 1/g^ = 0.75 at zero temperature, Eo — 0.346. A bit of arithmetic and the fact that 
the temperature is related to the temporal extent of the lattice A^,-, T = N^^, produces the phase diagram 
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of Figure 9. The solid line is the Mean Field phase boundary between the low-T, low-/i symmetry broken 
phase and the chiral symmetry restored phase. Note that the point at Tc/Sq = 0.482 and /ic/So = 0.564 
comes from the smallest lattice N-^ ~ 6 and is subject to the largest finite size correction. The data point at 
771 = comes from ref. [5] . 

It would be interesting to determine the order of the phase transition in Figure 9 away from the /i = 
and T = boundaries to test the Mean Field prediction of a line of second order transitions inside the plane. 
Our attempts to do this used time histories of S to search for two-state signals either visually or through 
histogramming. The results were not conclusive, however. In light of finite size effects to be discussed below, 
we believe that larger Nr values would actually be needed to accomplish this goal. 

6. Finite Size and Finite Mass Effects 

Finally, we did two additional simulations motivated as much by lattice QCD simulations as our interest 
in three dimensional four Fermi models. It is clear from the lattice Mean Field discussion in Section 3 above 
that even in the large Nf limit where fluctuations are suppressed relatively large lattices are needed to 
simulate the theory's actual critical behavior. We show this effect in Figure 10 which shows the induced 
fermion number plotted against ji for a 12"^ lattice at l/.g^ = 0.70. The discontinuous transition seen on a 
20'^ lattice in Figure 3 is replaced by a relatively smooth curve. A careful search for metastability in the 
region 0.35 < fi < 0.40 revealed none. 

It is also of interest to add a small explicit chiral symmetry breaking fermion bare mass to the Lagrangian 
and record its tendency to smooth out the transition. In fact, even a very small bare fermion mass obscures 
the first order transition of Figure 3 entirely. This effect is shown in Figure 11 where the theory is again 
simulated on a 20'^ lattice at l/.g^ = 0.70 and the induced fermion charge is measured with the bare mass 
of the dynamical fermion chosen to be either m ~ 0.01 or 0.005 in lattice units. A comparison with Figure 
3 shows that < n{ii,m) > decreases as m is increased from zero; in either phase this reflects the fact that 
more energy is always needed to excite fermion states from the vacuum once a bare mass is introduced. The 
average value of the sigma field is shown in Figure 12 for the same conditions. In neither case did we find 
any evidence for metastability, although, as shown in the figures, a very fine grid of fi values were simulated. 
One can also use the m = 0.01 and 0.005 results in Figure 12 to linearly extrapolate the average values of S 
to zero fermion mass at each value of /i. One can read off the figure that this procedure predicts a critical 
point at He = 0.42 — 0.43, considerably higher than the actual value of 0.3933(8) shown in Figure 1. A more 
sophisticated extrapolation method appears required to achieve quantitative results. QCD enthusiasts may 
find similar problems in their more complicated systems. 
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7. Conclusions 

In this study we have presented results of numerical simulations of an interacting relativistic field 
theory near its continuum limit at non-vanishing baryon-number density on far larger systems than have 
been possible hitherto [3] . The main result of our work is that the data we have is in excellent qualitative 
agreement with analytic predictions based on the 1/-/V/ expansion. We observe that in the regime \/ g^ < 1/gJ 
where chiral symmetry is spontaneously broken at zero temperature and density, for sufficiently large chemical 
potential fi there is a phase transition which restores the symmetry. For T = the transition is first order: for 
T > the transition appears continuous. Both results are consistent with the leading order l/Nf prediction 
that the transition at (/i/So = 1, T = 0) is an isolated first order point; however, the possibility of first order 
behavior persisting for values of T > cannot be excluded. In all cases our measured values of Tidi^T) fall 
slightly below the predictions of Sec. 2 - it will require further systematic study to establish whether this is 
due to finite volume effects or to genuine l/Nf corrections. 

It is interesting to contrast this situation with that of the two-dimensional Gross-Neveu model, where 
the leading order l/Nf expansion predicts a very rich phase diagram with a tricritical point in the (/i, T) 
plane separating first and second order symmetry restoring transitions [17]. Once quantum fluctuations are 
switched on, however, (ie. Nf < oo), the symmetry restoration is dominated by condensation of kink - anti- 
kink states, which are not present in the l/Nf expansion, and the transition becomes first order everywhere 
apart from an isolated second-order point at /i == [11]. 

We can speculate on whether the first order nature of the T = transition of the three-dimensional 
model remains stable as T is increased from zero, ie. whether there is in fact a tricritical point somewhere in 
the (/i, T) plane in this case, around which S exhibits power-law scaling rather than the essential singularity 
predicted in eqn. (3.13). Studies on much larger systems, enabling us to probe sufficiently low temperatures 
free of finite volume effects, would be required to locate such a point unambiguously. However, we note 
supporting evidence for the existence of such a point from a comparison on Figs. 3 and 10, showing plots 
of < n > vs. /i on 20"^ and 12"^ systems respectively. The major difference in the plots occurs in the broken 
phase, where < n > is considerably suppressed on the larger system, supporting the assertion made in Sec. 
4 that the signal in this phase is a finite volume effect (mean field theory predicts < n >= for T = in 
the broken phase, eqn. (3.14)). Even on the 20'^ lattice the presence of a signal shows that the system is 
experiencing a small but non-zero effective temperature: however the symmetry-restoring transition on this 
lattice is clearly first order, suggesting that the discontinuous nature of the transition is stable some way into 
the (/i,T) plane. Studies of the 0{l/Nf) corrections would be of value here: a first order transition would 
manifest itself in an extra non-trivial solution of the gap equation for fi,T > 0, implying an extra extremum 
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in the effective potential for S [17]. 

By studying the scaling of the critical chemical potential ^c as a function of coupling, we have been 
able to extract the correlation length critical exponent v. This is unusual, because we have only measured 
bulk properties of the system, and not any two-point correlations. Of course, our analysis implicitly relies on 
the theory's renormalizability [5,6]. Our result is in good agreement with the leading order 1/Nf prediction 
v — 1^ and our precision is slightly better than that obtained through finite-size scaling studies [5]. As always 
in these measurements, the major uncertainty is the location of the bulk critical point 1/ gl- Clearly from 
our data much greater precision will be required to probe 0{1/Nf) corrections. 

Finally, our studies on small volumes and with non-zero bare masses highlight the difficulties that must 
be faced when trying to understand the critical behavior of the chiral/thermodynamic limit by extrapolation 
from systems away from these limits. In this case the first order nature of the zero temperature transition 
is obscured - this should be a warning for the lattice QCD community (if one is still needed!) that a 
quantitative understanding of the behavior of QCD at non- vanishing density lies along a difficult road. 
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Table Captions 

1. Chemical Potential /i, ground state expectation value of the sigma field S, energy density e and ground 
state expectation value of the induced fcrniion number < n > on a 20'^ lattice at coupling 1/g^ = 0.70. The 
number of trajectories for the hybrid Monte Carlo algorithm is recorded in the last column. 

2. Same as Table 1 except l/g^ = 0.75. 

3. Measurement of S vs. ^ on Nr x N"^ [N = 36) lattices. The coupling l/g^ = 0.75. 
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Table 1 



M 


S 


£ 


<n> 


Trajectories 


.10 


.430(1 


.287(1) 


.00013 


2500 


.20 


.429(1 


.288(1) 


.00015 


2500 


.30 


.423(1 


.290(1) 


.0013 


2500 


.35 


.401(1 


.295(1) 


.0036 


1500 


.36 


.396(2 


.297(1) 


.005(1) 


1500 


.37 


.373(3 


.299(2) 


.006(1) 


1500 


.38 


.351(4 


.304(3) 


.008(1) 


1500 


.39 


.318(6 


.313(3) 


.013(1) 


3000 


.3925 


.307(8 


.316(3) 


.015(2) 


3000 


39375 


.275(8 


.321(3) 


.017(2) 


6000 


.3941 


.00 


.334(3) 


.026(2) 


6000 


.395 


.00 


.335(2) 


.026(2) 


3000 


.40 


.00 


.342(2) 


.029(2) 


1500 
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Table 2 



M 


S 


s 


<n> 


Trajectories 


.10 


.346(1) 


.300(1 


.00035(70) 


2500 


.20 


.340(2) 


.302(1 


.00032(65) 


1500 


.30 


.282(3) 


.312(2 


.006(2) 


1500 


.31 


.261(7) 


.316(2 


.009(1) 


1500 


.32 


.221(9) 


.322(2 


.012(1) 


3000 


3225 


.00 


.331(3 


.016(1) 


3000 


.325 


.00 


.330(3 


.015(1) 


3000 


.33 


.00 


.333(3 


.017(1) 


1500 


.335 


.00 


.336(3 


.020(1) 


3000 


.34 


.00 


.337(1 


.021(1) 


3000 



19 



Table 3 



iV^ = 6 Nr^8 Nr^lO Nr^U 






.259(2) 





.318(3) 





.335(2) 





.342(2) 


.10 


.236(3) 


.10 


.310(3) 


.20 


.307(2) 


.10 


.338(2) 


.15 


.181(3) 


.20 


.260(4) 


.25 


.258(3) 


.20 


.322(2) 


.16 


.167(3) 


.25 


.175(4) 


.26 


.243(3) 


.25 


.290(3) 


.17 


.140(3) 


.26 


.126(4) 


.27 


.202(3) 


.26 


.281(3) 


.18 


.090(8) 


.27 


.063(20) 


.28 


.192(3) 


.28 


.241(3) 


.19 


.050(15) 


.28 


.00 


.285 


.122(4) 


.30 


.138(7) 


.20 


.00 






.29 


.100(10) 


.31 


.030(15) 










.30 


.00 


.32 


.00 
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Figure Captions 

1. Mean Field predictions for E vs. fi curves for 1/g^ = 0.75. Solid lines are solutions for the contimuim 
Eq.(3.8) for Nr ~ 6, 8, 10 and 12 reading from the bottom left to the top right. The dotted lines are solutions 
of the discrete Eq.(3.18) (where convergent) for Nr = 6, 8, 10 and 12. The difference between the two sets 
of curves gives a measure of discretization effects. 

2. Chiral order parameters S vs. chemical potential /i on a 20'^ lattice at coupling l/g^ = 0.70. 

3. Induced ground state fermion number < n > vs. /x on a 20'^ lattice at l/g^ = 0.70. The dashed line is 
the mean field prediction. 

4. Chiral order parameters S vs. chemical potential /i on a 20'^ lattice at coupling l/g^ = 0.75. 

5. Induced ground state fermion number < n > vs. /i on a 20^ lattice at 1/g^ = 0.75. 

6. /ic/So vs. l/g"^ for a 20^ lattice. 

7. In^c vs. \n{l/gl - 1/g^) for l/g"^ = 0.70,0.75 and 0.80 on a 20^ lattice with the bulk transition 
1/gl = 0.975. 

8. S vs. /i for asymmetric lattices Nr x A^^ with A^^ ~ 6,8, 10 and 12 and A^ = 36 at coupling 1/g^ — 0.75. 

9. The phase diagram, ficf^a vs. Tc/Tiq, for the three dimensional Four-Fermi model. The solid Mean Field 
line separates the chirally broken phase at low ^ and T from the symmetric phase at large /i and T. 

10. The induced fermion number < n > vs. /i at coupling 1/g^ = 0.70 on a "small" 12^ lattice. 

11. < n > vs. ^ at 1/ g^ =■ 0.70 on a 20'^ lattice for the theory with a bare fermion mass term, ra = 0.01 and 
0.005 in lattice units. 

12. S vs. /i for the same parameters as in Fig. 11. 
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